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Abstract: We present an imaging method, dSLIM, that combines a novel 
deconvolution algorithm with spatial light interference microscopy (SLIM), 
to achieve 2.3x resolution enhancement with respect to the diffraction limit. 
By exploiting the sparsity of the phase images, which is prominent in many 
biological imaging applications, and modeling of the image formation 
via complex fields, the very fine structures can be recovered which were 
blurred by the optics. With experiments on SLIM images, we demonstrate 
that significant improvements in spatial resolution can be obtained by the 
proposed approach. Moreover, the resolution improvement leads to higher 
accuracy in monitoring dynamic activity over time. Experiments with 
primary brain cells, i.e. neurons and glial cells, reveal new subdiffraction 
structures and motions. This new information can be used for studying vesi- 
cle transport in neurons, which may shed light on dynamic cell functioning. 
Finally, the method is flexible to incorporate a wide range of image models 
for different applications and can be utilized for all imaging modalities 
acquiring complex field images. 
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1. Introduction 

Classical light microscopy techniques cannot be used directly in imaging most biological struc- 
tures, as they do not significantly absorb or scatter light [1]. Interference-based methods such 
as phase contrast [2] and differential interference contrast microscopy [3] allows imaging these 
transparent structures without the need for staining or tagging. Recently, more advanced meth- 
ods introduced the ability to measure quantitative information on the specimen by precisely 
quantifying optical phase shifts induced by the structure and motion of the specimen [4] . Spatial 
light interference microscopy (SLIM) [5], is a new and powerful quantitative imaging technique 
which allows high phase sensitivity imaging of nanoscale structures. SLIM has the important 
advantages of utilizing illumination with short-coherence length, and the ease of implementa- 
tion via add-on modules on existing phase-contrast microscopes. 

Although interference-based microscopy has tremendous advantages, it is still affected by 
the optical degradation and noise introduced by the instrument [6]. These degradations can be 
removed to a certain extent by employing post-processing methods. Deconvolution is a com- 
mon postprocessing method to invert the optical transfer function of the instrument. Although 
it is widely used in intensity-based microscopy [7-1 1], not much work has been reported on de- 
convolution in microscopy systems collecting quantitative information through complex fields. 
The work in [12] investigated the use of complex field deconvolution through inverse filtering 
in digital holographic microscopy [13], and have reported that the noise amplification, com- 
monly encountered with inverse filtering in intensity-imaging, is not as significant in the case 
of complex field microscopy. A nonlinear deconvolution method has been developed in [14] 
for SLIM that estimates the unknown magnitude and phase fields via a combination of variable 
projection and quadratic regularization on the phase component. 

In this paper, we present a novel method, dubbed dSLIM, for complex field deconvolution 
using an image model suitable for characterizing the fine scale structures. Based on the promi- 
nent features of phase images, we model the underlying image using the sparsity properties 
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of the transform coefficients. This model is especially useful in capturing fine-scale structures, 
and it successfully reveals the details in the phase components lost due to the instruments op- 
tical transfer function. In addition, due to the very low noise floor provided by SLIM (0.3 nm 
spatially and 0.03 nm temporally [5]), accurate experimental estimates of the point spread func- 
tion can be obtained, and deconvolution artifacts are significantly reduced. We demonstrate that 
resolution increases by a factor of 2.3 can be achieved with dSLIM. Thus, in dSLIM, images 
with a final resolution of 238nm can be rendered with only a 0.65NA objective. Moreover, the 
presented methodology is very flexible in incorporating image properties in a wide range of ap- 
plications (such as materials imaging) and can also be utilized for other interference microscopy 
techniques. 

The rest of this paper is organized as follows. We provide an overview of SLIM in Section 2. 
The image degradation model and the general framework for complex deconvolution is pre- 
sented in Section 3. The proposed deconvolution algorithm dSLIM is developed in Section 4. 
We present experiments with SLIM images in Section 5 and conclude in Section 6. 

We use the following notation throughout the paper: Bold letters h and H denote vectors and 
matrices, respectively, with transposes and H^. The spatial coordinates within a image are 
denoted by (x, j), operator * denotes convolution, and / is equal to \A-T. Finally, {•} is used to 
denote a set created with its argument. 

2. Overview of Spatial Light Interference Microscopy (SLIM) 

In SLIM [5], a spatially coherent light source, U {x^y) = \U{x^y)\ exp [/0(x,3;)], is used for illu- 
mination, which is decomposed into scattered and unscattered fields after passing the specimen. 
Let us denote the unscattered light as Uq and the scattered light as Ui{x^y). A liquid crystal 
phase modulator (LCPM) is used to introduce phase modulations to the unscattered field, such 
that 

f/o = |t/o|exp[-/(/)o], (1) 
Ui{x,y) = \Ui{x,y)\ cxp[i^i{x,y)] , (2) 

where (j)o is the intentionally added phase delay, and (j)i{x^y) is the phase difference between 
the scattered and unscattered fields caused by the specimen. The unscattered field contains the 
uniform background of the image field, whereas the scattered light provides information on the 
structure of the specimen. The recorded intensity is expressed as 

/(x,j,(/)o) = |t/oP + |t/i(x,3;)p + 2|[/o||[/i(x,3;)|cos(0i(x,j) + (^o). (3) 

In traditional phase-contrast microscopy [2], (^o is fixed at | and a single image is acquired. This 
can only provide qualitative information (i.e., (^i (x, cannot be uniquely retrieved). In contrast, 
SLIM uses multiple phase delays 0, |, ;r, and such that (j)i can be uniquely determined. 
Specifically, (j)i can be extracted from the four recordings using 



(pi{x^y) = arctan 



/(x,j,-f)-/(x,3;,f) 



I{x,y,0)-I{x,y,7r) 
Moreover, the phase associated with the complex field can be calculated by 

m{x,y) sm{(l)i {x,y)) 



0(x, j) = arctan 



l^m{x,y) sin((^i {x,y)) 



(4) 



(5) 



where we define by m{x^y) = the ratio of magnitudes of the scattered and unscattered 

light. 
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3. Image formation and deconvolution model 

As demonstrated in the previous section, both the magnitude and phase of the complex im- 
age function U{x^y) = j)| exp [/0(x, j)] can be uniquely determined using SLIM (using 
Eqs. (3) and (5)). However, as in all imaging systems, only a degraded version of this field can 
be observed in practice. Modeling the imaging process as a linear, spatially invariant degrada- 
tion system, the measured image can be expressed as the convolution of the original complex 
field with the instrument point spread function (PSF) as 

U{x,y) = U{x,y)^h{x,y)^n{x,y), (6) 

where h{x^y) is the PSF of the system, and n{x^y) is the additive signal independent noise. 
In general, both the magnitude and phase of the complex image function is degraded via the 
optical transfer function. As in traditional deconvolution [10], these fields can be estimated 
using a regularized inverse formulation 

\lJ{x,y)\M^,y)= argmin ^ \\U{x,y) -U{x,y)^h{x,y) f j)|,0(x, j)), 

(7) 

where is the noise variance, and the functional R{-) is used to regularize and impose con- 
straints on the estimates of the magnitude and phase. 

The estimation of both the magnitude and phase in Eq. (7) is a nonlinear optimization prob- 
lem, which is highly ill-posed and hard to solve in practice. However, in practice, the degrada- 
tion in the magnitude is very small compared to the degradation in the phase [12]. In addition, 
the phase image contains most of the information of the specimen relevant and useful from an 
application point of view. Therefore, it is convenient to assume that the magnitude of the field 
is constant, and the imaging system introduces negligible distortion in the magnitude, such that 
\U {x^y) \ ^ \U {x^y) \ = const. This assumption makes the problem in Eq. (7) linear, and it is also 
very useful in avoiding instabilities due to nonlinearity. With this approximation, the problem 
Eq. (7) becomes 

6(x,3;)= argmin || exp [/0(x,3;)] - /z(x,3;) * exp [/0(x,3;)] |p +j37?(0(x,3;)) . (8) 

For mathematical convenience and clarity, let us denote by g{x^y) the observed field 
exp [/<I)(x, j)] , and by f{x^y) the unknown field exp [/0(x, j)]. Due to the linearity of the degra- 
dation, the problem Eq. (8) can be expressed equivalently in matrix vector form as 

f=argmin^ ||g-Hf|p+j37?(f), (9) 

where g and f are images gix^y) and f{x^y) in vector forms, respectively, and H is the convo- 
lution matrix corresponding to the PSF h{x^y). 

When no regularization is used in Eq. (9), the closed- form solution of f can be found as 
(H^H) ^ H^g (equivalent to the inverse filter). However, this approach generally leads to noise 
amplification and ringing artifacts due to heavy suppression of high spatial frequencies. The role 
of regularization is to impose desired characteristics on the image estimates to avoid this noise 
amplification and to increase the resolution. The parameter j3 is used to control the trade-off 
between the data-fidelity and the smoothness of the estimates. 
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(a) (b) 



Fig. 1. Contourlets in (a) horizontal and (b) vertical directions. 



4. Complex field deconvolution using sparsity 

4.1. Image model 

It is well-known that phase-contrast imaging is highly sensitive to object boundaries but rela- 
tively insensitive to the flat background areas. Due to this characteristic, phase images exhibit 
high contrast around edges and local spatial variations within the specimen, which allows the 
capture of accurate shape and edge information. Moreover, in many important applications 
where interference microscopy is utilized, e.g., live cell imaging, the specimen contains very 
fine structures and small-scale movements that need to be precisely localized. 

Based on these observations, our goal is to construct a model of phase images that accurately 
represents this fine structure with sharp boundaries. We base our modeling on the sparsity 
principle, that is, our main assumption is that the phase images can be very accurately repre- 
sented in some transform domain with sparse coefficients [17, 18]. This transform sparsity can 
be achieved by appropriately selecting the transforms that capture the characteristics of spatial 
variations within the image. 

In this work, we consider a set of L linear transforms Dy^ of the complex image f with 
k = 1, . . . ,L. These transforms are chosen to be high-pass filters, such that their application 
provide complex images with a large number of coefficients with small values with only a few 
coefficients containing the most of the signal energy. The selection of the linear transforms 
that most accurately capture the image characteristics is crucial in the final image quality. We 
employ a collection of difference operators to capture signal variation at varying scales. The 
directional contourlets [19], depicted in Fig. 1, are used to capture the overall spatial variation, 
as they contain both vertical/horizontal and diagonal directions. In addition, for smaller scale 
features, we include first order directional difference operators 



and 45*^ and —45*^ first-order derivative filters 



" -1 0 " 




" 0 -1 " 


0 1 


7 


1 0 



(11) 



More complicated transforms can also be incorporated in the proposed framework in a straight- 
forward manner (possibly at the expense of computational complexity). 

Using these transforms, the image model can be constructed to exploit the sparsity in the 
transform coefficients. A commonly used sparse image model is [20] 



where || • ||p denotes the /^-pseudonorm, and are the weighting coefficients. It is known from 
the compressive sensing and sparse representation literature [15, 16, 20] that using 0 < /? < 1 
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enforces sparsity on Dy^f, with smaller p values increasing the sparsity effect. Hence, the prior in 
Eq. (12) enforces sparsity in the transform coefficients Dy^f, which in turn leads to smoothness 
in the image estimate. 

The disadvantage of using the prior in Eq. (12) is that it is nonconvex, and using this penalty 
creates a high number of local minima when estimating f. Instead, we use separate Gaussian 
priors on each transform coefficient 

p(f|{A,})ocexp(-itx«^ni(D^f)Jl2) , (13) 

\ k=\i=\ J 

or in a more compact form as 

p(f|{A,}) - exp [-\^ , (14) 

where are diagonal matrices with a^f, i = 1, . . .A/^ in the diagonal. Compared to Eq. (12), 
where a single parameter is assigned to all coefficients of filter output, separate parameters 
are used for each coefficient. It can be shown that Eq. (14) is equivalent to Eq. (12) in the limit 
p ^ 0 [21], hence it highly enforces sparsity. The model Eq. (14) has the advantage of being 
convex (as opposed to Eq. (12)), and therefore optimization over Eq. (14) is much easier and 
more robust compared to Ip minimization. 

The parameters a^i have a special important role in Eq. (13): they represent the local spatial 
activity at each location, and hence they are a measure of spatial variation in the correspond- 
ing filters direction. It is clear that the model Eq. (14) requires a large number of parameters, 
whose manual selection is not practical. We can, however, estimate them simultaneously with 
the complex image. For their estimation, we employ an additional level of model and assign 
uniform priors 

p{oCki) = const^ Mk^i. (15) 

Notice that this modeling assigns equal probability to all possible values of a^i, hence no prior 
knowledge is assumed on its value. 

It should be emphasized that this modeling based on sparsity principles does not necessarily 
cause the estimates to have very sparse coefficients. Real images are generally only approx- 
imately sparse, i.e., a few transform coefficients have large values whereas most coefficients 
are very small, but not necessarily exactly zero. This behavior is generally referred to as com- 
pressible [22]. Enforcing sparsity to an extreme extent can therefore suppress subtle image 
features, which may be important. Our modeling in Eq. (14), on the other hand, allows for the 
small-valued transform coefficients through Gaussian distributions while enforcing the general 
compressible structure of the images. 

4.2. Noise model 

The signal-independent noise is modeled via a independent Gaussian noise model on the ob- 
served field as 

p(g|f,a2)ocexp(^-^||g-Hf||2^, (16) 

with the noise variance. An additional level of modeling (as in Eq. (15)) can be incorporated 
to estimate this parameter as well. However, SLIM provides images with very high SNRs (on 
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the order of 1000 or more), and therefore O"^ is generally very small and can be estimated 
experimentally from an uniform area of the observed image. In addition, the Gaussian noise 
assumption becomes an accurate description of the noise in SLIM due to the high SNR (as in 
fluorescence microscopy [10]). 

4.3. Algorithm 

Using the modeling described in the previous sections, we formulate the problem of estimating 
the unknown complex image f and the parameters a^i using the maximum a posteriori (MAP) 
estimates, that is. 



f , d;t, = argmin -log 



p(gif,a2)p(f|{A,})nnpK) 

k=\i=\ 



(17) 



1 ^ 

= argmin || g - Hf ||^ + ^ (^k^)" {^k^) • (18) 

f,«fc.- ^ k=\ 

This problem is convex in f and a^i, but not jointly, and therefore we resort to an iterative 
scheme to estimate the unknowns in an alternating fashion. The optimal estimate of the complex 
image can be found by taking the derivative of Eq. (18) and setting it equal to zero, which results 
in 

f = I^H^H + X D^A^D^ j H^g . (19) 

The parameters a^i can be estimated in a similar way by equating the corresponding deriva- 
tives to zero, which results in 

^ki = -—^2 ' ^20) 

where e is a small number (e.g., 10~^) used to avoid numerical instability. It follows from 
Eq. (20) that the parameters a^i are functions of the filter response at pixel /, and therefore 
a spatially -adaptive estimation is employed for f in Eq. (19) through their joint estimation. 
Notice also that matrices Ay^ are spatial-adaptivity matrices controlling the smoothness applied 
at each location; when the filter responses at a pixel are very small, the algorithm assumes that 
the pixel has low spatial variation, and applies a large amount of smoothness at that point. On 
the other hand, if the filter responses are high, the pixel is likely to be close to an edge and the 
smoothness amount is lowered to preserve the image structure. 

In summary, the proposed method dSLIM consists of alternating estimations of the complex 
field f using Eq. (19), and the spatial adaptivity matrices Ay^ using Eq. (20). The block diagram 
of a single dSLIM iteration is shown in Fig. 2. The filters Dy^ consist of the derivative operators 
in Eq. (10) and Eq. (11), and the directional contourlets shown in Fig. 1. The estimate of the 
image Eq. (19) can be computed very efficiently using the conjugate gradient (CG) method. The 
matrices H and T>k do not have to be explicitly constructed during CG iterations; all operations 
in Eq. (19) can be performed via convolutions in the spatial domain or multiplications in the 
Fourier domain. Empirically, we found that the algorithm converges rapidly; a few iterations 
(up to 5-10) is generally enough to provide high-quality results. Hence, the proposed method 
can be applied to large images very efficiently. 

The proposed algorithm contains only one free parameter (the noise variance G^) that needs 
to be set by the user. In our experiments, we empirically estimated its value by taking a rect- 
angular region of the image with uniform values and computing the variance in this region. As 
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Fig. 2. Block diagram of one dSLIM iteration. 

mentioned above, this estimate is known to be reliable in images with high SNR [10]. Since 
this free parameter corresponds to a physical quantity, its estimation is relatively easy and does 
not require extensive image-dependent tuning. 

Finally, it should be emphasized that dSLIM does not alter the quantitative imaging property 
of SLIM, which is one its main advantages. As the deconvolution is applied to the complex field 
exp [/<I)(x, j)] rather than to the measured intensities /(x, (j)o) in Eq. (3), the quantitative phase 
information is preserved. In contrast, traditional deconvolution methods [10] applied directly to 
the intensity images can not preserve the quantitative information. 

5. Experiments 

In this section, we illustrate the application of dSLIM to complex field images obtained by 
SLIM and quantitatively demonstrate the resolution increase. All SLIM images were acquired 
using a white-light source (mean wavelength A = 530 nm); the field of view is 75/im x 100/im 
with the CCD resolution of 1040 x 1388. In all reported experiments, the specimen is relatively 
thin such that the whole image is in focus, and the degradation in the image is only due to 
a planar PSF. The PSF, depicted in Fig. 3(a), is obtained experimentally by imaging a sub- 
resolution 200nm microbead treated as a point- source. Due to the high SNR provided by SLIM, 
this PSF closely matches the actual optical transfer function of the imager. 

In all images, the noise level is estimated within the range 10~^-10~^ (for a maximum signal 
value of 1), which is used as the value of the parameter O"^. The NAs of the objective and 
condenser are NAq = 0.75 and NAc = 0.55, respectively. The experimentally measured full- 
width-at-half-maximum (FWHM) of the PSF is 540 nm, which is comparable with the expected 
Rayleigh limit, calculated as j^^'^j^^ = 497nm. 

We first investigate the resolution increase obtained by dSLIM by applying it to the exper- 
imental PSF. The experimental PSF is shown in Fig. 3(a). Treating this as the original image, 
we apply dSLIM and obtain the result shown in Fig. 3(b). The FWHM of the original PSF is 
approximately 540nm, whereas after dSLIM, the FWHM is reduced to approximately 238nm, 
corresponding to a 2.3 times increase in resolution. The horizontal cross-section of the images 
are shown in Fig. 3(c), where the reduction in FWHM is clearly visible. Notice that this res- 
olution is significantly below the diffraction limit. The estimated FWHM during the iterative 
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Fig. 3. (a) Experimental PSF, (b) result of dSLIM, (c) normalized horizontal cross-sections, 
and (d) estimated FWHMs during the iterative procedure. 
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Fig. 4. Images of two microbeads (a) SLIM and (b) dSLIM. The cross-sectional profiles 
(with normalized maximum phase values) are shown in (c). 



procedure is shown in Fig. 3(d), which shows a significant reduction in the first iterations and 
convergence within 10 iterations. 

To further examine the resolution increase provided by dSLIM, we next apply it to SLIM 
images of multiple 200nm microbeads. Figure 4(a) shows an image of two beads approximately 
550 nm apart, which are barely resolved in the original SLIM image. The image after applying 
dSLIM is shown in Fig. 4(b), where the microbeads are clearly separated while their distance 
is accurately preserved (Fig. 4(c)). 

Next, we demonstrate dSLIM images of biological specimen. A SLIM phase image of a 
hippocampal neuron is shown in Fig. 5. The SLIM images are shown on the left column, while 
dSLIM images are shown on the right column. It is clear that dSLIM effectively removes the 
blur, while deconvolution artifacts and noise are successfully suppressed and object boundaries 
faithfully preserved. dSLIM recovers the details of the fine structure of the specimen which 
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Fig. 5. (a) SLIM image of a hippocampal neuron, (b) image provided by dSLIM. Detailed 
images of the central parts are shown in (c) and (d). 



are hard to observe in the original image. The high quality of the dSLIM result can be better 
observed from the detailed parts of the image shown in Fig. 5(c) in comparison with the SLIM 
image. 

Another example is shown in Fig. 6 (Media 1, Media 2), which is a snapshot of a dynamic 
SLIM image sequence of a live hippocampal neuron culture. Interference-microscopy is ex- 
tremely useful in monitoring dynamic cellular processes over time, as it does not require in- 
vasive contrast enhancement techniques (such as fluorescence tagging). SLIM is a very attrac- 
tive modality for this application due to its very high spatial and temporal resolution (several 
frames/second). It is clear from Fig. 6(c) (left) (Media 2) that the observation noise level is very 
low, but the image exhibits a certain level of blur due to the diffraction-limited PSF. The result 
of dSLIM is shown in Fig. 6(c) (right) (Media 2), which shows a clear resolution improvement 
over the original image. The increased spatial resolution also positively affects the examina- 
tion of dynamic neuron processes. Due to more accurate estimates of size and locations of the 
particles, the dynamic changes and hence the biological behavior can be better observed (see 
Media 2 for a visualization). 

To further confirm the increase in resolution in real images, we examine microparticles in the 
SLIM image shown in Fig. 6. Figures 7(a) and 7(b) show detailed images of a single particle 
from the SLIM and dSLIM images (marked as region D in Fig. 6(a)) . The cross-sections of 
the images are shown in Fig. 7(c). The vertical diameter of the particle is measured as approx- 
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Fig. 6. SLIM dynamic imaging of live hippocampal neuron in primary cell culture, (a) 
SLIM image (Media 1), (b) dSLIM image (Media 1), (c) Detailed areas of regions A, B, 
C from (a) and (b) (Media 2). The SLIM image regions are shown on the left, while the 
dSLIM image regions are on the right. 

imately 1.5/im in the SLIM image, whereas it is measured as 0.63/im in the dSLIM image. 
The reduction in the length of the particle is approximately 2.3, which is in agreement with the 
result of the PSF deconvolution experiment (Fig. 3). 

Our final experiment shows two neuronal processes (putative axons) which were not resolved 
in the original SLIM image (region marked as E in Fig. 6(a)). The detailed area is shown 
in Fig. 8(a), and the dSLIM result is shown in Fig. 8(b). dSLIM reveals two objects located 
approximately 430 nm apart. This can also be observed from the normalized cross- sections 
passing through the maximum phase values, shown in Figs. 8(c). The evolution of this area over 
time is shown in Fig. 9 (Media 3). It can be observed both from the original images and the 
cross-sections that the objects are just resolved in some time frames, but unresolved in others. 
dSLIM successfully separates the objects through the whole dynamic sequence (Media 3). 
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Fig. 8. (a) Two very closely located particles from the region E in Fig. 6 not resolved 
in the SLIM image, (b) dSLIM image, and (c) the normalized cross-sections through the 
maximum phase values. 





Fig. 9. Dynamic evolution of the area shown in Fig. 8 (Media 3). Top row: SLIM images, 
middle row: dSLIM images, bottom row: Normalized cross-sections of the images (through 
the segment shown in the top-left image) at each time point. 



#144841 - $15.00 USD Received 28 Mar 201 1; revised 4 May 201 1; accepted 24 May 201 1; published 2 Jun 201 1 
(C) 201 1 OSA 1 July 201 1 / Vol. 2, No. 7 / BIOMEDICAL OPTICS EXPRESS 1826 



6. Conclusion 



In this paper, we presented a novel deconvolution method, dSLIM, for complex image fields 
acquired by interference microscopy. Our formulation is based on three key observations: First, 
the image formation can be treated as a linear process in the complex fields, such that the degra- 
dation of the microscopy can be modeled by a PSF acting on the complex images. Second, due 
to the high SNR provided by the SLIM, the experimentally obtained PSF of the imager closely 
matches the actual PSF. Finally, the phase images of biological specimen can be very accurately 
modeled using sparsity principles. We combined these properties to develop a very effective de- 
convolution procedure that significantly improves the final resolution, allowing imaging very 
fine structures and motions in live cells below the diffraction limit. Due to the high spatial and 
temporal resolution, this approach can be utilized to acquire new information for studying live 
cells. 
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